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Abstract 

In this paper we develop further a method for detecting unstable pe- 
riodic orbits (UPOs) by stabilising transformations, where the strategy is 
• to transform the system of interest in such a way that the orbits become 

^ ^ stable. The main difficulty of using this method is that the number of 

r— H transformations, which were used in the past, becomes overwhelming as 

, ^ , we move to higher dimensions 1191 120| . We have recently proposed a set 

of stabilising transformations which is constructed from a small set of al- 
^— ( ready found UPOs 2 . The main benefit of using the proposed set is that 

^ its cardinality depends on the dimension of the unstable manifold at the 

UPO rather than the dimension of the system. In a typical situation the 
dimension of the unstable manifold is much smaller than the dimension 
of the system so the number of transformations is much smaller. Here we 
extend this approach to high-dimensional systems of ODEs and apply it to 
7—i the model example of a chaotic spatially extended system - the Kuramoto- 

T-H Sivashinsky equation. A comparison is made between the performance of 

this new method against the competing methods of Newton- Armijo (NA) 
and Levernberg-Marquardt (LM). In the latter case, we take advantage of 
^ the fact that the LM algorithm is able to solve under-determined systems 

of equations, thus eliminating the need for any additional constraints. 

X 

1 Introduction 

The following work is concerned with the detection of UPOs for large systems 
of ODEs of the form 



with X e M" and v{x) G M". Often Eq. (1.1) is the result of a spatial dis 



cretisation of a parabolic PDE. In this way dynamical system results can be 
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applied to extended systems which exhibit chaos. Over the past twenty years 
the periodic orbit theory (POT) has been developed and successfully applied to 
low dimensional systems. Many dynamical invariants such as natural measure, 
Lyapunov exponents, fractal dimensions and entropies |18j can be determined 
via periodic orbit (a.k.a cycle) expansions. It is an open question whether or not 
the POT has anything to say for spatially extended systems. It is thus an im- 
portant numerical task to find UPOs for such systems in an attempt to answer 
this question. The model example of a spatiotemporally chaotic system is the 
Kuramoto-Sivashinsky equation (KSE) , which was first studied in the context of 
reaction-diffusion equations by Kuramoto and Tsuzuki |12j . whilst Sivashinsky 
derived it independently as a model for thermal instabilities in laminar flame 
fronts [13] . It is one of the simplest interesting PDEs to exhibit chaos and we 
have chosen it as a test case for the work that follows. 

The problem of finding UPOs is essentially a root finding problem, thus a 
popular approach is to use some variant of Newton's method. Indeed, Zoldi and 
Greenside have reported the detection of 127 distinct UPOs for the KSE 
Here they solve for the discretised system using simple shooting with a damped 
Newton method to update each step. Another recent assault on the KSE by 
Cvitanovic and Lan use variational methods [3l [13] . Here a suitable cost func- 
tion is constructed so that its minimisation leads to the detection of a UPO. 
The main problem, however, is that Newton type methods suffer from two ma- 
jor drawbacks: firstly, with increasing period, the basins of attraction become so 
small that placing an initial seed within the basin is practically impossible, and 
secondly, the method has no way of differentiating between true roots and local 
minima of the cost function. The latter drawback is one which increases sig- 
nificantly with dimension due to the complicated topology of multi-dimensional 
fiows. 

The method of detecting UPOs by stabilising transformations [Hl[51[Tl] aims 
at transforming the system in such a way that its UPOs become stable. Unlike 
in the Newton-type methods, the transformations are linear and thus do not 
suffer from spurious convergence. When faced with the task of finding UPOs of 
a discrete system 

x,+i^f{x,), /:M"^M", (1.2) 
one can look instead at the related fiow 

^=9i^), (1-3) 
as 

where g{x) = P{x) — x. It is straightforward to see that the period— p points 
of the map are equilibrium points for the associated flow. With this setup we 



are able to stabilise all UPOs x* of Eq. (1.2 1 such that all the eigenvalues of 
the Jacobian DfP{x*) have real part smaller then one. In order to stabilise all 
possible UPOs we study the following flow 
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where C S K"^" is a constant matrix introduced in order to stabilise UPOs with 
the Jacobian's that have eigenvalues with real parts greater than one. Given a 
set {C} of such matrices, we have a family of differential equations which need to 
be solved in order to find all UPOs of Eq. ( 1.2 1. One example of such a set was 
proposed by Schmelcher and Diakonos (SD) I19J. It is the set Csd of orthogonal 
matrices such that only one entry {±1} per row or column is nonzero. It has 
been verified that the set Csd stabilise all hyperbolic fixed points for n < 2, and 
numerical evidence suggests the result holds for n > 2, but, thus far, no proof 
has been presented. However, if we wish to extend the stabilising transformation 
approach to higher dimensions, the set Csd cannot be applied directly, since its 
size increases very rapidly with system dimension (|Csd| = 2"n!). 

For high-dimensional systems with relatively few unstable directions, the 
method of stabilising transformations can be applied efficiently by restricting 
our attention to the unstable part of Eq. (1.2 1. Indeed, by constructing trans- 
formations which only alter the stability of the flow in Eq. ( 1.3 1 in the unstable 
subspace of DfP, it is possible to reduce the number of transformations consid- 
erably. The authors have recently proposed a new set of matrices C based on 
the properties of a small already detected set of UPOs. Here the cardinality is 
|C| = 2"", where n„ is the dimension of the unstable manifold at x* . This is 
the key to extending these ideas to higher dimensional systems, since often in 
practice the systems of interest are such that n„ <C N. For example, after a 
finite difference discretisation of the KSE with resulting system of size N — 100, 
only four of the corresponding Lyapunov exponents are positive [35]. Thus at 
each seed we would have only |C| = 16 matrices as opposed to |Csd| = 2^^991 if 
we used the SD matrices 



2 Subspace decomposition 



In what follows we take our leave from the subspace iteration methods [151 E2] 
Consider the solution of the nonlinear system 



0, 



/ 



(2.1) 



where /(x) is assumed twice differentiablc in the neighbourhood of x* , an iso- 
lated root of Eq. (2.1 ). We can approximate the solution of (2.1 1 by a recursive 



fixed point procedure of the form 



Xi+i^fixi), i= 1,2,3,. 



(2.2) 



It is well known that the iteration (2.2 1 converges locally in the neighbourhood 
of a solution x* , as long as all the eigenvalues of the Jacobian Df{x*) lie within 
the unit disc {z G C : |z| < 1}. In contrast, (2.2 1 typically diverges if Df{x*) 



has an eigenvalue outside the unit disc. In that case, a popular alternative is to 



^Here we reduce the flow to a Poincare surface of section in order to obtain a discrete 
system 
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employ Newton's method 



{Df{x,)-I^)5x, = (2.3) 
Xi-\-i — Xi -\- dXi-j i — 1,2,3,.... 

The idea of subspace iterations is to exploit the fact that the divergence of 
the fixed point iteration (2.2 1 is due to a small number of eigenvalues, n„, lying 
outside the unit disc. By decomposing the space M" into the direct sum of the 
unstable subspace spanned by the eigenvectors of Df{x*) 

P = Span{efe G M" : Df{x*)ek - XkCk, \Xk\ > 1} 

and its orthogonal complement, Q, a modified iterative scheme is obtained. The 
application of Newton's method to the subspace P whilst continuing to use the 
relatively cheap fixed point iteration on the subspace Q, results in a highly 
efficient scheme provided dim(P) ^ dim(Q). 

To this end, let Vp E M"^"" be a basis for the subspace P C M" spanned 
by the eigenvectors of Df(x*) corresponding to those eigenvalues lying outside 
the unit disc, and Vq € M"^"" a basis for Q, where n„ + = n. Then, we can 
define orthogonal projectors P and Q onto the respective subspaces, P, Q, as 
follows 



Note that any x G 
X = Vpp - 



Q = VgVj = I,,-P. (2.4) 
" admits the following unique decomposition 
Vqq^P+q, p:^Vpp^Px, q --Vqq^ Qx, (2.5) 



with p e M"" and q e M"=. Substituting (2.5 1 in Eq. (2.3) and multiplying the 
result by [Vq, Vp]^ on the left, one obtains 



Vq^DfVq - In 







DfVq Vp' DfVp - /„„ 



Aq 

Ap 



■ q 
P 



(2.6) 



' p ' 1 — '-'riu xris I y q " p — "rig xn^ ; and DfVp 



Here we have used the fact that VpVq = On,,xns i ^J^p = ^r, 
On,xn„ the latter holding due to the invariance of Df on the subspace P. Now, 
the first equations in (2.6 1 may be solved using the following fixed point 
iteration scheme 



Aq = 



0, 

VjDfVqAq^' 



+ VJf - q, 



(2.7) 



/-I 



Aq^'^=y2{VjDfVqy{Vjf-q), 



=Q 



where I denotes the number of fixed point iterations taken per NR step. Since 
rc[Vj DfVq] < 1 by construction, the iteration (2.7 1 will be locally convergent 
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on Q in some neigbourhood of Aq 
order to determine Ap one solves 



here 



denotes the spectral radius. In 



{VjDfVp - /„ J Ap = -Vjf + p- VjDfVgAq. 



Note that in practice only one iteration of (2.7) is performed [TF, i.e. I = 1, this 
leads to the following simplified system to solve for the correction [Aq, Ap]^ 



-In, 

VjDfVg VjDfVp 



Aq 
Ap 



vJf 



Key to the success of the above algorithm is the accurate approximation of 
the eigenspace corresponding to the unstable modes. In order to construct the 
projectors P, Q, the Schur decomposition is used. However, primary concern of 
the work in [151 122] is the continuation of branches of periodic orbits, where it 
is assumed that a reasonable approximation to a UPO is known. Since we have 
no knowledge apriori of the orbits whereabouts we shall need to accommodate 
this into our extension of the method to detecting UPOs. 



2.1 Stabilising transformations 

As discussed in the introduction, an alternative approach is supplied by the 



method of STs, where in order to detect equilibrium solutions of Eq. (2.1 1 we 
introduce the associated flow 



dx 
ds 



= Cg{x). 



(2.8) 



Here g = fP{x) — x an d C G 



is a constant matrix. 



Now, substituting (2.5 1 in Eq. (2.8) and multiplying the result by [V^,!^]^ 
on the left, one obtains 



dq 
ds 
dp 

ds 



V^g. 



(2.9) 
(2.10) 



Thus we have replaced the original Eq. (2.8) by a pair of coupled equations 



Eq. (2.9) of dimension Us and Eq. (2.10) of dimension n„. Since 

and r„ \Vj D gV^^] < by construction, it follows that in order to detect all UPOs 
of Eq. (12. 11) , it is sufficient to solve 



dq 

ds 
dp 

ds 



= v: 



cv^g. 



(2.11) 
(2.12) 
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(a) Schur decomposition 
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X 

(b) Singular value decomposition 



Figure 1: Basins of attraction for the period— 3 orbits of the Ikeda map with 
parameter values a — 1.0, b = 0.9, k = 0.4 and rj = 6.0. Here we have chosen 
C = — 1 since in this example the unstable subspace is one-dimensional. 



where C G M""^*"" is a constant matrix. 

In [151 [22] the Schur decomposition is used in order to construct the projec- 
tors P and Q. This is fine for continuation problems since one may assume from 
the offset that they possess an initial condition xo sufficiently close to a UPO 
such that the Schur decomposition of Df{xo) gives a good approximation to the 
eigenspace of Df{x*). However, it is well known that the eigenvectors of the 
perturbed Jacobian Df(x* + 5x) behave erratically as we increase Sx. In order 
to enlarge the basins of attraction for the UPOs we propose that singular value 
decomposition (SVD) be used instead. That is we choose an initial condition xq 
and construct the SVD of its pre-image, i.e. Df{f^^{xQ)) = USW''^ (or in the 
continuous case Defy'" (xq)) — USW''^ for some time T), the columns of U 
give the stretching directions of the map at xq , whilst the singular values deter- 
mine whether the directions are expanding or contracting. It is these directions 
which we use to construct the projectors P and Q. Due to the robustness of 
the SVD we expect the respective basins of attraction to increase. 

in order to apply the 



It is not necessary in practice to decompose Eq. (2 



new ST. Rather we can express C in terms of C and Vp. To see this we add Vq 
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times Eq. (|2.11| to Vp times Eq. (|2.12[) to get 



= VgV;^gix) + VpCV;^g{x), 
= [In + Vp{C-I^jV;]g{x), 



dx 
ds 



where the second hne follows from Eq. ( |2.4[ ) . From this we see that the following 
choice of C is equivalent to the preceding decomposition 

C = /„ + ^p(C-/„jF/. (2.13) 

Thus in practice we compute Vp and C at the seed xq in order to construct C 
and then proceed to solve Eq. ( 2 



The advantage of using the SVD rather than the Schur decomposition can 
be illustrated by the following example. Consider the Ikeda map [H]: 



/(x) 



a + b{x^ cos (0j) - sin (4)) 
b{xi sin + j/i cos {(j)^)) 



(2.14) 



where (pi = k — ri/(l + xf + yf) and the parameters are chosen such that the map 
has a chaotic attractor: a — 1.0, b — 0.9, k — 0.4 and rj — 6.0. For this choice of 
parameters the Ikeda map possesses eight period— 3 orbit points (two period-3 
orbits and two fixed points, one of which is on the attractor basin boundary). In 
our experiments we have covered the attractor for Eq. (2.14 1 with a grid of initial 
seeds and solved the associated flow for p = 3, i.e., g{x) — f^{x) — x. This is 
done twice, firstly in the case where the projections P and Q are constructed via 
the Schur decomposition and secondly when they are constructed through the 
SVD. Since all UPOs of the Ikeda map are of saddle type, the unstable subspace 
is one-dimensional and we need only two transformations: C — I and C = — 1. 
Figure [l] shows the respective basins of attraction for the two experiments with 
C = — 1. It can be clearly seen that the use of SVD corresponds to a significant 
increase in basin size compared to the Schur decomposition. Note that with 
C = —1 we stabilise four out of eight fixed points of f^. The other four are 
stabilised with C = 1. The corresponding basins of attraction are shown in 
Figure |2] Note that for this choice of ST the decision of which basis vectors to 
choose is not important, since Eq. (2.13) yields C = /, so that the associated 



flow is independent of the selected basis. 



3 Implementation 

A typical approach in the determination of UPOs for flows is via a Poincare 
surface of section (PSS). By "clever" placement of an (n— l)-dimensional man- 
ifold in the phase space, the problem is reduced to a discrete map defined via 
intersections with the manifold. However, a correct choice of PSS is a chal- 
lenging problem in itself. Due to the complex topology of a high-dimensional 
phase space, the successful detection of UPOs will be highly dependent upon 



7 



-1 



1 
X 



2 



Figure 2: The basins of attraction for the Ikcda map for the choice of C = 
1. Fixed points of with negative unstable eigenvalues are stable stationary 
solutions of the associated flow, while those with positive eigenvalues are saddles 
located at the basin boundaries. 



the choice of surface. When the choice of a suitable PSS is not obvious a priori, 
we found it preferable to work with the full flow, adding an auxiliary equation 
to determine the integration time T. 

Let X 1-^ (f>*{x) denote the flow map of Eq. ( 1.1 1. Then we define the associ- 
ated flow as follows 

^ = C{cl,^{x)-x). (3.1) 
as 

The additional equation for T is constructed such that T is always changing in 
the direction that decreases \(j)^{x) — x\, i.e. 



ds 



d\9l 
dT ' 



or, more precisely 



dT 
ds 



-av{<j>^{x)) ■ {(l)'^{x) - x). 



(3.2) 



Here a > is a constant which controls the relative speed of convergence of 
Eq. (3.2 1. This leads to the following augmented flow which must be solved to 
detect UPOs of (O): 



d 


X 




ds 


T 





C{(j)^{x) - x) 
-av{(j)^ (x)) ■ {(t)^{x) — x) 



(3.3) 



Note that the augmented flow (3.3) is derived by integrating a nonlinear PDE 
for some time T and will become increasingly stiff for larger T. 

Several approaches have been proposed for the solution of stiff systems of 
ODEs; see for example, the review by Shampine and Gear [2T]. Of all these 
techniques the general-purpose codes contained within the ODEPACK software 
package [7] are regarded as some of the best available routines for the solution 
of such systems. Thus, in our numerical experiments we use the stiff solver 



dlsodar from the ODEPACK toolbox to integrate (3.3 1, dlsodar is a variable 



step-size solver which automatically changes between stiff and nonstiff modes. 



In particular, as we approach a steady state of Eq. (3.3 1 dlsodar will take 
increasingly larger time-steps, leading to super linear convergence in the neigh- 
bourhood of the solution. 

To use the solve r dl sodar, we must provide a routine that returns the value 
of the vector field ( |3.3[ ) evaluated at a given point (a;,T)^. Here we need the 
solution of Eq. ( |1.1[ ) which is obtained by applying a suitable numerical integra- 
tion scheme; see the next section for further details. The ODEPACK software 
package makes use of the Jacobian matrix of the system being integrated and 
provides the option of computing the Jacobian via finite differences or via a 
user supplied routine. Note that in the case that the flow is expected to be stiff 
much of the time, it is recommended that a routine for the Jacobian is supplied 



and we do this. The derivative of (3.3) with respect to {x,T)~^ is given by 



C( Jt - /„) CvT 
— a(uj( Jt — In) + 9^DvJt) —a{g^ Dvvt + v^vt) 

where Jt = dcf)^ {x) /dx, vt = Dv — dv/dx and g = cj)^{x) — x a,s 

usual. 

Quite often one might wish to terminate simulation before the usual stopping 
criteria of, for example, a maximum number of steps taken or certain tolerances 
having been reached. A nice feature of the dlsodar algorithm is that it al- 
lows the optional user supplied routine to do just this. To be more exact, it 
determines the roots of any of a set of user supplied functions 

hi^hi{t,Xi,...,Xn), i=l,...,m. 



and returns the solution of (3.3 1 at the root if it occurs prior to the normal 
stopping criteria. 

Note that, to increase the efficiency of the algorithm we wish to avoid the 
following two instances: flrstly, due to the local nature of the STs, we should 
stop the search if we wander too far from the initial condition, and secondly. 



since our search is governed by the dynamics of Eq. (3.1 1 and not by those of 



( 1 . 1 1 , we might move off the attractor after some time period so the convergence 
to a UPO becomes highly unlikely. In our numerical experiments we supply the 
following function 

hi = a - \\g\l 

where a £ M is a constant and || • || denotes the L2 norm. In practice, we have 
found that there exists a threshold value of a, such that convergence is highly 
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unlikely once the norm of g surpasses it . Note that we also restrict the maximum 
number of allowed integration steps since the convergence becomes less likely if 
the associated flow is integrated for a long time. 

We must also provide two tolerances, rtol and atol, which control the local 
error of the ODE solver. In particular, the estimated local error in X = {x, T)^ 
will be controlled so as to be less than 

rtol • \ \X\\oo + atol. 

Thus the local error test passes if, in each component, cither the absolute error 
is less than atol or the relative error is less than rtol. The accuracy with which 



we would like to solve the flow (3.3) is determined by the stability properties 
of To understan d th is, we note that in evaluating the RHS of (3.3 1 it 



is the solution of Eq. (1.1 1 at time T, i.e. (j>^{x), which is critical for error 
considerations. Suppose that our initial point lies within 6x of a true trajectory 
X. Then (/)^(a; + Sx) lies approximately within e^'^Sx of the true trajectory, 
(j)^{x), where A is the largest Lyapunov exponent of the system. Since A is 
positive for chaotic systems, the error grows exponentially with the period and 
we should take this into account when setting the tolerances rtol and atol. This 
leads us to the following settings for the tolerances 

rtol = atol = 10" Ve^"^". 

where Tq is the initial period and A is the largest Lyapunov exponent of the 
flow V. We have computed the Lyapunov exponent using the algorithm due to 
Benettin ei aZ [1]. 

3.1 Kuramoto-Sivashinsky equation 

We have chosen the Kuramoto-Sivashinsky equation (KSE) for our numerical 
experiments. It is the simplest example of spatiotemporal chaos and has been 
studied in a similar context in [U [T31 US], where the detection of many UPOs 
has been reported. We work with the KSE in the form 

ut^^-{u^).-u,.-u,..., (3.4) 

where x E [0,i] is the spatial coordinate, t E M+ is the time and the sub- 
scripts X, t denote differentiation with respect to space and time. For L < 2tt, 
u{x, t) = is the global attractor for the system and the resulting long time dy- 
namics are trivial. However, for increasing L the system undergoes a sequence 
of bifurcations leading to complicated dynamics; see for example [TT]. 

Our setup will be close to that found in [I^]. In what follows we assume 
periodic boundary conditions: u{x,t) — u(x + L,t), and restrict our search to 
the subspace of antisymmetric solutions, i.e. u(x,t) — —u{L — x,t). Due to 



the periodicity of the solution, we can solve Eq. (3.4 1 using the pseudo-spectral 



method [HJIH]. Representing the function u{x,t) in terms of its Fourier modes: 

u{x,t) ■.^T-^[u]^Y.Uke-'^'i\ 

fcez 
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where 



1 

, UO,Ul,...)^, Uk J^[u]k = -jr J U{x,t)t 



we arrive at the following system of ODEs 



duk 
~dt 



ikq . 



Here q — 271 / L is the basic wave number. Since u is real, the Fourier modes are 



related by u_ 



Furthermore, since we restrict our search to the subspace 



of odd solutions, the Fourier modes are pure imaginary, i.e. £He(ufe) = 0. 

The above system is truncated as follows: the Fourier transform T is replaced 
by its discrete equivalent 



7V-1 



Ofc := TN\u]k 



J=0 



N-l 
k=0 



where Xj = L/N and aif^k = o,*k- Since — Q due to Galilean invariance and 
setting ajv/2 = (assuming N is even), the number of independent variables in 
the truncated system is n = N/2 — 1. The truncated system looks as follows: 



dk = [{kqY - {kqY]ak 



ikq 



(3.5) 



with k — 1, . . . , n, although in the Fourier transform we need to use au over the 
full range of k values from to — 1. 

The discrete Fourier transform Tn can be computed using fast Fourier trans- 
form (FFT). In Fortran and C, the routine REALFT from Numerical Recipes [25] 
can be used. In Matlab, it is more convenient to use complex variables for a^. 
Note that Matlab function f f t is, in fact, the inverse Fourier transform. 

To derive the equation for the matrix of variations, we use the fact that Tn 
is a linear operator to obtain 



dhk 

dan 



[{kqf-{kqY]5k,+kqTN[T],\a]®T^\5k,]], J = 1, . . . , TV - 2, (3.6) 



where ® indicates componentwise product, and the inverse Fourier transform is 
applied separately to each column of 5kj- Here, 5kj is not a standard Kronecker 
delta, but the N x n matrix: 



/ 
1 








V-1 
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with index k running from to A'^ ~ 1. 

In practice, the number of degrees of freedom n should be sufficiently large 
so that no modes important to the dynamics are truncated, whilst on the other 
hand, an increase in n corresponds to an increase in computation. To determine 
the order of the truncation in our numerical experiments, we initially chose n to 
be large and integrated a random initial seed onto the attractor. By monitoring 
the magnitude of the harmonics an integer k was determined such that aj < 10~^ 
for j > k. The value of n was then chosen to be the smallest integer such that: 
(i) n > k, and (ii) N = 2n + 2 was an integer power of two. The second 
condition ensures that the FFT is applied to vectors of size which is a power of 
two resulting in optimal performance. 

Note that in the numerical results to follow we work entirely in Fourier space. 
We use an exponential time differencing method (ETDRK4) due to Kassam 



and Trefethen [9] in order to solve (3.51 and (3.6 1. Note that the method uses 



a fixed step-size (h = 0.25 in our calculations) thus it is necessary to use an 
interpolation scheme in order to integrate up to arbitrary times. In our work 
we have implemented cubic interpolation |25j . More precisely, to integrate up 
to time t e [ti,ti + h], where the ti are integer multiples of the step-size h. 
We construct the unique third order polynomial passing through the two points 
a(ti) and a{ti + h), with derivatives a'(ti) and a'(ti + h) at the respective points. 
In this way we obtain the following cubic model: 

p{s) = [2a{U) + a'{t,) + a {U + h) - 2a{U + h)]s^ + [3a{U + h) - 3a{U) 
- 2a'{U) - a'{ti + h)]s^ + a'{U)s + a{ti), 

where the parameter s — {t — ti) / h E [0, 1]. 
3.2 Numerical Results 

We now present the results of our numerical experiments. We have compared 
the performance of our method against the Newton- Armijo (NA) version of the 
damped Newton algorithm |10j, as well as the nonlinear least squares solver 
Imder from the MINPACK software package [17]. Note that Imder is an im- 
plementation of the Levenburg-Marquardt algorithm [12]. Both methods have 
been successfully applied to spatially extended systems in order to detect pe- 
riodic orbits in the past, in particular, in [26 the NA algorithm was able to 
detect many distinct UPOs of the KSE, whilst more recently, Imder has been 
used to determine many UPOs of the closely related complex Ginzburg-Landau 
equation [11] . 

Both methods require as input a function whose zeros are to be determined. 
In the case of NA an auxiliary equation must be added due to the time invariance 



of Eq. (1.1). For Imder, it is not necessary to augment the system since the 
algorithm is able to solve under-determined systems of equations, however, we 
run two separate experiments: (i) we solve the unconstrained system, and (ii) 
we supply an auxiliary equation as in the case of the Newton algorithm. In 
this way we are able to infer the effect of the PSS on the performance of the 
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search. One important note concerning Imder, is that, even though it is capable 
of solving under-determined systems of equations, it still requires the number 
of equations to equal the number of unknowns in the user supplied function, in 
practice however, we may simply set any additional equations identically equal 
to zero. 

In order to determine UPOs via the aforementioned methods, we introduce 
the following augmented system 



F{x,T) 



(t>^{x) — X 
v{xo) ■ {(t)^{x) - xq) 



(3.7) 



where the additional equation defines a Poincare surface of section normal to 
the initial field vector. Also, in order to optimise efhciency we use the analytic 
Jacobian in all our experiments rather than a numerical approximation 



DF{x, T) = 



Jt — In 



vt 




(3.8) 



Note that, in the case of the unconstrained system, the equations are as above. 



except that we set the last entry in (3.7) and the final row in (3.8 1 identically 



equal to zero. In addition to this we must supply the two routines with certain 
tolerances in order to control errors. Let us denote the relative error desired in 
the approximate solution by xtol and the relative error desired in the sum of 
squares by ftol. In the computations performed in the next section we set the 
tolerances to the recommended values of 



ftol = 10" 



xtol = 10" 



Further, Imder requires a third tolerance, gtol, which measures the orthogonality 
desired between the vector function F and the columns of the Jacobian, DF; we 
also set this to its recommended value gtol= 0. Finally, we specify a maximum 
number of function evaluations allowed during each run of the Imder and NA 
algorithms in order to increase efficiency. 

For our method we use the set of matrices proposed in §2.1| with 



C — CsD, 

since within the low-dimensional unstable subspace it is possible to apply the 
full set of Schmelcher-Diakonos (SD) matrices. The UPOs determined from our 
search will then be used as seeds to determine new cycles. Here we proceed 
in analogous fashion to P| by constructing STs from the monodro my m atrix, 
DffF (a*), of the cycle {a* ,T*Y . We then solve the augmented flow (3.3 1 from 

,r)^, where the time T is chosen such that 



the new initial condition (a* 

a*(0) w a*(f) 



and f (modT*) ^ 0. 



Note that any given cycle may exhibit many close returns, particularly longer 
cycles, thus in general a periodic orbit may produce many new initial seeds. 
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This is an especially useful feature, since we do not have to recompute the 
corresponding STs. 

Initially, to determine that a newly detected cycle, (a;*,r*), was different 
from those already found, we first checked whether the periods differed, that is, 
we determined whether \T* —T'\ > Ttoi for all previously detected orbits. If two 
orbits where found to have the same period, we then calculated the distance 
between the first components of x*, and all other detected orbits y* . However, 
in practice we found, that if two orbits have the same period, then either they 
are the same or they are related via symmetry, recall that if u{x, t) is a solution 
then so is —u{L — x,t). Thus, in order to avoid the convergence of UPOs that 
are trivially related by symmetry we will consider two orbits as being equal if 
their periods differ by less than the tolerance Ttoi- Of course, this criterion can, 
in theory, lead us to discard cycles incorrectly, however, this is highly unlikely 
in practice. 



3.2.1 Comparison of the numerical methods 

The search for UPOs is conducted within a rectangular region containing the 
chaotic invariant set. Initial seeds are obtained by integrating a random point 
within the region for some transient time r. Once on the attractor, the search 
for close returns within chaotic dynamics is implemented. That is, we integrate 
the system from the initial point on the attractor until a(io) ~ o,{ti) for some 
times to < ti, and use the close return, (a(io),'7o), where Tq = ti — to, as our 
initial guess to a time-periodic solution. In order to build the STs we solve the 
variational equations for each seed starting from the random initial point, a(0), 
for time r + to- In order to construct the matrix Vp, we apply the SVD to the 
matrix 

D(t>^+'°{a{0)) ^USW^, (3.9) 
the corresponding ST is given by 

where Vp — Ujk, j — 1, . . . , n, k — 1, . . . , n^, i.e. the first columns of U in 



(3.9 1. Here is the number of expanding directions which is determined by 
the number of singular values with modulus greater than one. 

We examine two different system sizes: L = 38.5, for which the detected 
UPOs typically have one positive Lyapunov exponent, and L = 51.4, for which 
the UPOs have either one or two positive Lyapunov exponents. The correspond- 
ing systems sizes are n = 15 and rt = 31 respectively. Our experiments where 
conducted over two separate ranges. We began by looking for shorter cycles 
with period T € [10, 100], the lower bound here was determined a posteriori so 
as to be smaller than the shortest detected cycle. We then searched for longer 
cycles, T S [100, 250] to be more precise, where the maximum of T = 250 was 
chosen in order to reduce the computational effort. 



In our calculations we set the positive constant a = 0.25 in Eq. (3.3 1. Using 



the solver dlsodar we integrated 500 random seeds over both ranges for time 
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Table 1: The number of distinct periodic solutions for the Kuranioto-Sivashinsky 
equation detected by the method of STs. Here L — 38.5 and a = 0.25. 



Period 


C 










Work 




Co 


28 


498 


252 


10 


412 


10-100 




16 


296 


684 


32 


1196 






44 


794 


936 


42 


1608 




Co 


235 


395 


903 


78 


2151 


100 - 250 




64 


256 


1294 


112 


3086 






299 


651 


2197 


190 


5237 



Table 2: The number of distinct periodic solutions for the Kuramoto-Sivashinsky 
equation detected by the Levernberg-Marquardt algorithm Imder with L — 38.5. 



Method 


Period 






A^fcv 




Work 


Constrained 


10 - 100 
100 - 250 


41 
232 


475 
358 


12 
113 


9 

97 


156 
1665 


Unconstrained 


10 - 100 
100 - 250 


42 
233 


491 
363 


11 
109 


9 

95 


155 
1629 



s = 150, the seeds where chosen such that \(j)^°{a{0)) — a(0)| < 1.0. If the 
flow did not converge within 1000 integration steps, we found it more efficient 
to terminate the solver and to re-start with a different ST or a new seed. As 
mentioned in the previous section we choose a constant a = 50 experimentally 
so that integration is terminated if the norm of g grows to large, i.e. |(7(a;)| = 
\(j)^{x) — x\ > a. Typically on convergence of the associated flow the UPO 
is determined with accuracy of about 10^*" at which point we implement two 
or three iterations of the Newton- Armijo rule to Eq. (3.7 1 in order to allow 
convergence to a UPO to within roundoff error. 

Similarly, we run the Imder and NA routines from the same 500 seeds over 
the two different time ranges. The respective routines are terminated if one of 
the following three scenarios arise: (i) a predefined maximum number of function 
evaluations is exceeded, we set the maximum number of function evaluations 
equal to 1000, (ii) the error between two consecutive steps is less than xtol, but 
the sum of squares is greater than ftol, indicating a local minimum has been 
detected, or (iii) both xtol and ftol are satisfied indicating that convergence to 
a UPO has been obtained. 

Note that, one problem with applying methods that use a cost function to 
obtain "global" convergence to our setting, is that for increasing period, the 
level curves of \g\'^ become increasingly compressed along the unstable manifold 
of (fF {x), resulting in a complicated surface with many minima, both local and 
global, embedded within long, winding, narrow "troughs" . 
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Table 3: The number of distinct periodic solutions for the Kuramoto-Sivashinsky 
equation detected by the Newton- Armijo algorithm with L = 38.5. 



Period 


A^po 








Work 


10-100 


43 


475 


29 


7 


141 


100 - 250 


213 


334 


309 


41 


965 



This can be explained by the following heuristics: for simplicity let us assume 
that we are dealing with a map Xn+i = f{xn)i whose unstable manifold is a 
one-dimensional object. In that case, we may define a one-dimensional map, 
locally, about a period-p orbit, a;*, of the map / as 

h{s) = \g{x* + 6x)\\ 

Here g = fP{x) — x as usual, 6x = x{s) — x* is small, and we only allow x{s) 
to vary along the unstable manifold. Expanding g in a Taylor series about the 
periodic orbit, x*, we obtain 

h{s) = \g{x*) + Dg{x*)6x + 0{Sx^)\'^ 
= \{DfP{x*)^I^)Sx + 0{Sx^)\^ 
= 6x'^{AP - InfSx + 0{Sx^), 

where the third line follows since 6x is an eigenvector of Df{x*) with correspond- 
ing eigenvalue A, and A = diag(A, ... A). Note that in the above we assume that 
p is large but finite so that the term (A^* — remains bounded. Now, close to 
the periodic orbit, h is approximately a quadratic form with slope of the order 
A^, and since |A| > 1, it follows that if we move along the unstable manifold 
I^P will grow quicker for larger periods, or, in other words, the level curves are 
compressed along the unstable direction. 

Since both Imder and NA reduce the norm oi g = 4>'^{x) — a; at each step, 
they will typically follow the gradient to the bottom of the nearest trough, where 
they will start to move along the narrow base towards a minimum. Once at the 
base of a trough, however, both algorithms are forced into taking very small 
steps, this follows due to the nature of the troughs, i.e. the base is extremely 
narrow and winding, and since both methods choose there next step-size based 
on the straight line search. In order to avoid this situation we use the following 
additional stopping criteria in our experiments 



-^Vstep > 



MXn)\\ 



(3.10) 



Here iVstop denotes the maximum number of iterations allowed, and the term 
on the right is a linear approximation of the number of steps required for con- 
vergence. If the above condition fails in 50 consecutive iterates we terminate 
the search. 
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Figure 3: Illustration of two UPOs of KSE detected from a single seed. We 
show both a level plot for the solutions and a projection onto the first two 
Fourier components. Since u{x,t) is antisymmetric on [0,L], it is sufficient to 
display the space-time evolution of u{x,t) on the interval [0,L/2]: (a) Seed 
with time T = 37.0, (b) a periodic solution of length T = 36.9266 detected 
with stabilising transformation C = +1 and (c) a periodic solution of length 
T — 25.8489 detected with stabilising transformation C — —1. 



In order to make a comparison between the efficiency of the methods we 
introduce the measure of the work done per seed 

Work = iVfev + nx iVjov. (3.11) 

Here A^fcv is the average number of function evaluations per seed, Nj^v the 
average number of Jacobian evaluations per seed and n is the size of the system 



being solved. The expression in (3.11 1 takes into account the fact that evaluation 
of the Jacobian is n times more expensive than evaluation of the function itself. 

The results of our experiments are summarised in Tables [l]"|6] Here A'po 
denotes the number of distinct orbits found, A'hit gives the number of times we 
converged to a UFO, and Nf^v, N^cv, and Work, are as defined above. In Tables 
[T] and |4] the performance of the stabilising transformations are analysed both 
collectively and on an individual basis; here the different Ci denote the different 
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Table 4: The number of distinct periodic solutions for the Kuramoto-Sivashinsky 
equation detected by the method of stabilising transformations. Here L — 51.4 
and a = 0.25. 



Period 




AT 

JVpo 


AT 

iVhit 


AT 

JVfcv 


AT 


Work 




Co 


11 


365 


382 


16 


894 




Ci 


1 


111 


411) 


23 


114d 






7 


lUo 


A 


9n 
zu 


luyo 




'-'3 


9 


oi 


00 / 


10 


77'X 
t 1 


10-100 


C4 


2 


157 


464 


26 


1296 




C5 





171 


654 


35 


1774 




Ce 


5 


174 


666 


36 


1818 






2 


139 


496 


36 


1648 






30 


1417 


3885 


205 


10445 




Co 


51 


330 


628 


17 


1172 




Ci 


7 


209 


807 


27 


1671 




C2 


17 


138 


936 


31 


1928 




C3 


21 


177 


917 


26 


1749 


100 - 250 


Ci 


11 


161 


877 


36 


2029 




C5 





116 


960 


43 


2336 




Ce 


1 


117 


975 


42 


2319 




C7 


6 


161 


879 


35 


1999 






114 


1409 


6979 


257 


15203 



SD matrices. In Table [T] Co = +1 and Ci = —1, whilst in Table [4] we have 



Ca 



1 

1 

1 

1 



Ci = 



,C5 



-1 
1 

-1 ■ 

1 



,^2 = 



, Cg 



1 

-1 

1 ■ 

-1 



,^3 = 



,C7 



-1 








-1 





-1 


-1 





The total work done per seed is given by the sum over all the SD matrices and 
is denoted by {Ci}. 

In total we found 487 UPOs using the ST method, 379 UPOs using Imder 
and 350 UPOs using NA. Whilst all methods found roughly the same number 
of orbits when searching for shorter period cycles. Tables [T]- [4] reveal significant 
differences in performance. By comparing the work done per seed we see that 
both Imder and NA are considerably faster, indeed efficiency between ST and 
the other methods for L = 38.5 differ by factors as great as 10. The situation 
changes, however, when we look at the detection of longer cycles. For L = 38.5 
in particular, we see that the ST method computes many more orbits than its 
competitors. Also, although both Imder and NA remain more efficient than the 
ST method, the difference in efficiency is now only a factor of 3 in the case of 
Imder and 5 in the case of NA. In fact, using the identity matrix alone, the 
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Table 5: The number of distinct periodic solutions for the Kuranioto-Sivashinsky 
equation detected by the Levernberg-Marquardt algorithm Imder with L = 51.4 



Method 


Period 






A^fev 




Work 


Constrained 


10-100 


19 


311 


26 


18 


602 


100 - 250 


85 


337 


72 


61 


2024 


Unconstrained 


10 - 100 


19 


489 


16 


12 


400 


100 - 250 


85 


396 


63 


51 


1695 



Table 6: The number of distinct periodic solutions for the Kuramoto-Sivashinsky 
equation detected by the Newton- Armijo algorithm with L = 51.4. 



Period 


^po 




A^fcv 




Work 


10 - 100 


19 


308 


73 


20 


713 


100 - 250 


75 


309 


174 


43 


1550 



ST method detects more UPOs than either Imder or NA, but with comparable 
efficiency. For larger system size, L = 51.4, the ST method still detects more 
orbits than Imder or NA. However, in doing so a considerable amount of extra 
work is done. It is important to note here, that the increase in work is due, 
mainly, to the fact that not all the SD matrices work well. For example, the 
subset of matrices {Cq, C2, C3, C4} C Csd, detect approximately 90% of the 
longer period UPOs, yet they account for less than 50% of the overall work. 
Note that it is not surprising that the SD matrices do not perform equally well, 
this is, after all, what we would have expected based upon our experience with 
maps. However, it is still important to try all SD matrices - when possible - in 
order to compare their efficiency, especially if one wishes to construct minimal 
sets of ST matrices. 

Another important consideration is, how was the performance of Imder af- 
fected by the additional constraint? We can see from Tables |3] and [s] that the 
unconstrained system is more efficient, in all cases requiring fewer steps to con- 
verge, more importantly, it can be seen that the constrained system can fail to 
converge at all. Crucially, the number of searches which failed increases consider- 
ably as we look at larger system size, therefore, as we move to more complicated 
systems we would expect to see this difference in performance further increase. 
For example, in [14J Lopez et al use Imder to search for relative periodic orbits 
in the complex Ginzburg-Landau equation, where they have augmented the sys- 
tem with three additional equations. Our results suggest that this search would 
have benefited, not only in performance and numbers of UPOs detected, but by 
the savings in both time and effort required to construct additional equations 
and the resulting Jacobian, by setting all additional equations identically equal 
to zero. 
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Finally, a key feature of our method is that we can converge to several 
different UPOs from just one seed, depending upon which ST is used. Figure |3] 
shows one of such cases, where Eq. (3.3 1 is solved from the same seed for each 
of the 2"" STs (ri„ = 1 in this example) , with two of them converging to two 
different UPOs. Figure [3^ shows the level plot of the initial condition and a 
projection onto the first two Fourier components. Figures [3]d and [3]: show two 
unstable spatiotemporally periodic solutions which where detected from this 
initial condition, the first of period T = 36.9266 was detected using C = +1, 
whilst the second of period T = 25.8489 was detected using C = —1. The ability 
to detect several orbits from one seed increases the efficiency of the algorithm. 



3.2.2 Seeding with UPOs 

In order to construct a seed from an already detected orbit, (a*, T*), we begin 
by searching for close returns. That is, starting from a point on the orbit, a* (to), 
we search for a time ti such that a{ti) is close to a{to). As long as TmodT 7^ 0, 
where T — ti — to, we take (a*(io),T) as our initial guess to a time periodic 
solution. For longer period cycles, we can find many close returns by integrating 
just once over the period of the orbit. Shorter cycles, however, produce fewer 
recurrences and, in general, must be integrated over longer times to find good 
initial seeds. In our experiments we searched for close returns T £ (10, 2T*) if 
T* < 150.0; otherwise, we chose f G (10, T*). 

Stabilising transformations are constructed by applying the polar decompo- 
sition to the matrix G = QB, where G is defined by 

G := y(S'A-I„)F-i. 

Here V and A are defined through the eigen decomposition of D(f>^ {a* (to)) = 
VAV~^, and S = diag(±l, ±1, . . . , ±1). The different transformations are then 
given by C = —Q^. Now, in ^ it was shown that a change in the signs 
of the stable eigenvalues does not result in a substantially different stabilising 
transformation, thus we use the subset of S such that Sa — 1 for i > 71^. For 
L — 38.5, we have just two such transformations, since all periodic orbits have 
only one unstable eigenvalue. Whilst for L = 51.4, we will have either two 
or four transformations depending upon the number of unstable eigenvalues of 
D</>^'(a*(to)). 

In the calculations to follow, a close return was accepted if \a*{ti) — a*{tf))\ < 
2.5. Note that, this was the smallest value to produce a sufficient number of 
recurrences to initiate the search. A particular cycle may exhibit many close 
returns. In the following, we set the maximum number of seeds per orbit equal 
to 5. To obtain convergence within tolerance 10~^ we had to increase the 
integration time to s = 250 and the maximum number of integration steps 
allowed to 2000. If we converged to a UPO then as in the preceding section, we 
apply two or three Newton-Armijo steps to the Eq. ( |3.7| in order to converge 
to machine precision. 

The results of our experiments are given in Tables [7] ancjS] As in the previous 
section, A^po denotes the number of distinct orbits found, A^hit the number of 
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Table 7: Number of distinct orbits detected using the method of stabihsing 
transformations with periodic orbits as seeds The number of seeds is 489. L = 
38.5, a = 0.25. 



c 




iVhit 


A^fcv 




Work 


Co 


46 


209 


2814 


136 


4854 




52 


198 


2819 


130 


4769 




98 


407 


5633 


266 


9623 



Table 8: Number of distinct orbits detected using the method of stabihsing 
transformations with periodic orbits as seeds. The number of seeds is 123. 
L ^ 51.4, a = 0.25. 



c 


Npo 


NMt 


^fcv 




Work 


Co 


1 


21 


2797 


48 


4333 


Ci 


2 


37 


2815 


40 


4095 










216 


3 


312 







1 


187 


2 


251 




3 


59 


6015 


93 


8691 



times we converged to a UPO, A^fev the average number of function evaluations 
per seed, and iVjac the average number of jacobian evaluations per seed. The 
computational cost per seed is measured in terms of the average number of 
function evaluations per seed and is defined as in Eq. ( [3lT| ). In Tables [7] and |8] 
the different Ci can be uniquely identified by the signature of pluses and minuses 
defined through the corresponding matrix S. For L = 38.5, the matrices Cq 
and Ci correspond to the signatures (+, ...,+) and (—,+,...,+), respectively, 
whilst for L = 51.4, the matrices Co, Ci, C2 and C3 correspond to (+, . . . , +), 
(— , +,..., +), (— , —,+,..., +), and (+, —,...,+) respectively. As before, the 
total work done per seed is defined as the sum over all matrices and is denoted 
by {C;}. 

For L = 38.5 we where able to construct 489 seeds from 343 previously 
detected periodic orbits. From the new seeds we found a further 98 distinct 
UPOs, bringing the total number of distinct orbits detected for L — 38.5 to 
433. An important observation to make is that both matrices perform equally 
well, as can be seen from Table [T] This is in contrast to the SD matrices for 
which the performance of the different matrices varies greatly; see Tables [T] and 
|4j The result of using periodic orbits as seeds has, however, increased the cost 
per seed as compared with the cost of seeding with close returns within a chaotic 
orbit, see Tables [l] and [Tj The increased computation is due mainly to the fact 
that the close returns obtained from periodic orbits are not as good as those 
obtained from a chaotic orbit. Recall that the stabilising transformations are 
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Figure 4: Illustration of how a UPO can be used as a seed to detect new cycles. 
We show the projection onto the first two Fourier components of the initial seed 
(T = 237.4470) and the three detected orbits. Here L = 38.5 and n = 15. 



based on the local invariant directions of the orbit, and we would expect the 
performance to suffer as we move further from the seed. 

For larger system size, where the system becomes more chaotic, the construc- 
tion of seeds from close returns becomes increasingly difficult. For L = 51.4, we 
constructed 123 initial seeds from the 144 periodic orbits using the method of 
near recurrences. Here we detected only 3 new distinct UPOs. It is important 
to note, however, that although we do not find many new UPOs for the case 
L — 51.4, the method is still converging for approximately 20% of all initial 
seeds. Also, in Table |8] the poor performance of the matrices C2, C3 as com- 
pared to that of Co, Ci, is due to the fact that only a small percentage of seeds 
were constructed from orbits with two positive Lyapunov exponents. 

One advantage of using periodic orbits as seeds is that we can construct many 
seeds from a single orbit. Figure |4] gives an illustration of three periodic orbits 
which were detected from a long periodic orbit used as seed. Figure |4] shows the 
projection onto the first two Fourier components of the initial seed (of period 
T = 237.4470) and the three orbits detected whose periods in descending order 
are T = 200.428, T = 58.7515 and T = 52.1671. For each seed determined from 
a particular periodic orbit the stability transformations are the same. Hence, 
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the ability to construct several seeds from a single orbit, increases the efficiency 
of the scheme. 

4 Conclusions and further work 

We have presented a scheme for detecting UPOs in high dimensional chaotic 
systems based upon the stabilising transformations proposed in [2, 5 , Due 
to the fact that one often wishes to study low dimensional dynamics embed- 
ded in a high dimensional phase space, it is possible to increase the efficiency 
of the stabilising transformations approach by restricting the construction of 
such transformations only to the low-dimensional unstable subspace. Following 
the approach often adopted in subspace iteration methods [T5|, we construct 
a decomposition of the tangent space into unstable and stable orthogonal sub- 
spaces. We find that the use of SVD to approximate the appropriate subspaces 
is preferable to that of Schur decomposition, which is usually employed within 
the subspace iteration approach. As illustrated with the example of the Ikeda 
map, the decomposition based on SVD is less susceptible to variations in the 
properties of the tangent space away from a seed and thus produce larger basins 
of attraction for stabilised periodic orbits. Within the low-dimensional unstable 
subspace, the number of useful stabilising transformations is relatively small, so 
it is possible to apply the full set of SD matrices. In fact, we have found that 
the subset of diagonal matrices of ±1 is capable of locating a large number of 
UPOs, although more analysis will be carried out in the future to determine if 
it is possible to detect all types of UPOs with this subset. 

In conclusion, we have presented an extension of the stabilising transfor- 
mations approach for locating periodic orbits in high-dimensional dynamical 
systems. Future work will concentrate on rigorous mathematical analysis of 
this approach in order to determine the range of its applicability. We will also 
work on the development, within the stabilising transformations approach, of an 
efficient strategy for systematic detection of periodic orbits in high-dimensional 
systems. 
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